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INTRODUCTION 


1 . 


This  Is  the  final  report  of  a  three  year  effort  to  study  physical 
processes  of  relevance  to  the  mass  spectrometric  measurement  of 
stratospheric  ions.  The  effort  involves  the  development  of  a  Monte  Carlo 
■odel  of  the  free  jet  expansion  occurring  within  the  mass  spectrometer 
including  the  effects  of  aggloaeration  onto,  and  fragmentation  of,  ionic 
clusters.  This  model  has  been  named  the  EXPANDO  code. 

The  attempt  to  carry  out  in  situ  mass  spectrometry  in  the  stratosphere 
is  complicated  by  collision  induced  changes  which  may  occur  in  the  gas 
stream  as  it  expands  after  passage  through  the  orifice.  Both  positive  and 
negative  ions  exist  In  the  stratosphere  with  clustered  polar  molecules 
surrounding  the  ion  core.  As  these  ion  clusters  are  carried  along  in  the 
expanding  gas  stream,  the  falling  temperature  will  tend  to  favor  the 
formation  of  larger  clusters.  The  charge-dipole  interaction  is 
characterized  by  large  cross  se(tions,  so  agglomeration  of  polar  molecules 
can  potentially  alter  thp  cluster  size  distribution  tlr  t  the  quadrupole 
sees  from  the  distribution  that  exists  in  the  undisturbed  stratosphere. 
Conversely,  the  measured  cluster  size  distribution  may  be  driven  towards 
smaller  clusters  via  fragmentat  on.  As  the  ionic  clusters  are  selectively 
accelerated  by  the  electric  fie  d  within  the  mass  spectrometer,  high  energy 
collisions  with  neutrals  may  bn  ak  apart  the  clusters. 

The  present  effort  invo  ves  a  Monte  Carlo  simulation  of  these 
processes,  so  that  a  model  can  be  used  to  relate  the  measured  properties  to 
those  existing  in  the  undisturbed  atmosphere.  The  direct  simulation  Monte 
Carlo  method  involves  storing  a  discrete  number  of  molecules  (via  their 
velocities,  positions  and  other  pertinent  information)  in  a  computer.  The 
solution  region  is  broken  up  into  a  number  of  cells,  and  the  solution  is 
stepped  forward  in  time  in  a  two  stage  process.  First,  the  molecules  are 
advanced  along  their  trajectories  by  an  amount  appropriate  to  their 


velocity  and  a  tiae  increment  &tB.  In  this  i irst  stage  soae  molecules  will 
leave  the  solution  region  and  soae  will  be  i  itroduced  as  determined  by  the 
boundary  conditions.  The  second  stage  is  to  simulate  collisions  in  each 
cell  appropriate  to  Atm  so  that  collision  frequencies  are  properly 
simulated.  A  basic  hypothesis  of  the  method  is  that  if  the  tiae  step  is 
made  small  enough  the  processes  of  translations  and  collisions  can  be 
uncoupled  in  this  manner. 

After  a  steady  state  has  been  achieved,  the  solution  is  periodically 
sampled  by  accumulating  statistical  sums  of  number  densities,  velocities 
and  other  basic  properties.  The  solution  is  run  repeatedly  until 
statistical  deviations  are  reduced  to  a  desired  limit,  and  then  physically 
meaningful  output  quantities  are  computed  from  the  statistical  sums.  The 
number  of  molecules  represented  is  typically  a  few  thousand  at  a  time, 
which  is  vastly  fewer  than  the  number  occurring  in  virtually  all  real 
flows.  Hence,  the  construction  of  a  dynamically  similar  flow  to  be 
simulated  in  the  computer  is  an  essential  feature  of  the  method. 

The  logic  of  the  solution  procedure  for  steady  state  flows  is  shown  in 
Figure  1,  which  includes  the  steps  described  above.  Section  2  is  a 
discussion  of  the  gas  dynamics  of  a  freejet  expansion,  and  Section  3 
presents  the  calculational  grid  developed  for  the  EXPANDO  code.  The  basic 
elements  of  the  direct  simulation  Monte  Carlo  method  have  been  described 
elsewhere.1,2  and  they  are  presented  here  only  briefly.  The  application  of 
the  technique  to  the  mass  spectrometer  problem  is  discussed  in  Section  4, 
and  the  chemical  scheme  for  proton  hydrate  reactions  is  developed  in 
Section  5.  Sample  calculations  and  results  are  presented  in  Section  6. 

1Bird,  G.  A.,  Molecular  Gas  Dynamics,  Clarendon  Press,  Oxford,  1976. 

p 

cElgin,  J.  B.,  "Monte  Carlo  Calculations  of  Mass  Spectrometer  Flow",  Report 
AFGL  TR-83-0057,  Air  Force  Geophysics  Labor, itory,  February  1983. 


Figure  1 


A  Diagram  of  the  Basic  Solution  Proi edure  utilized  in  the 
EXFAKDO  Monte  Carlo  Transition  Flow  Code. 


THE  FREEJET  EXPANSION  WITHIN  A  MASS  SPECTROMETER 


2  . 

2.1  Motivation  for  Studying  Major  Species  Freejet 

The  freejet  expansion  of  the  Major  (neutral)  species  through  a  sonic 
orifice  into  a  near  vacuus  is  a  classic  problem  representing  the  zeroth 
order  solution  to  the  problem  at  hand.  The  numerical  dominance  of  the 
major  species  assures  that  their  distribution  in  phase  space  (i.e.,  their 
density  and  velocity  distributions)  will  be  negligibly  affected  by  the 
■inor  species  of  interest  (e.g..  ion  clusters,  H20,  etc.).  Hence,  for  a 
given  atmospheric  pressure  and  orifice  geometry  the  freejet  expansion  of 
the  major  species  can  be  calculated  once  and  for  all.  The  physical 
processes  of  interest  (ion  acceleration  due  to  electric  fields, 
agglomeration,  fragmentation,  etc.)  are  then  handled  by  considering  the 
minor  species  to  be  traveling  within  a  known  phase  space  distribution  of 
major  species.  Considerable  conceptual  and  computational  simplifications 
result  from  the  recognition  that  the  freejet  expansion  is  uncoupled  from 
the  ionic  motion. 

2.2  Physical  Description  of  the  Freejet  Expansion 

The  major  features  of  a  freejet  created  by  expanding  air  through  a 
sonic  orifice  into  a  near  vacuum  are  illustrated  in  Figure  2.  For  the 
cases  of  interest  the  orifice  is  always  considerably  larger  than  an  ambient 
mean  free  path,  so  the  initial  portion  of  the  expansion  is  best  described 
in  terms  of  continuum  fluid  mechanics.  The  orifice  forms  a  sonic  throat  in 
which  the  flow  is  accelerated  to  a  Mach  number  of  unity.  If  the  orifice 
has  a  finite  thickness  ("t”  in  Figure  2),  then  a  boundary  layer  forms  on 
the  edge  of  the  orifice,  restricting  the  J  low  somewhat  from  that  which 
would  be  predicted  via  one-dimensional  in\iscid  theory.  The  mass  flow 
through  the  orifice  can  be  represented  by^ 
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where  rQ  is  the  orifice  radius.  T  is  the  ratio  of  specific  heats  in  the  gas 
and  P0  and  p0  are  the  stagnation  pressure  and  density,  respectively.  The 
discharge  coefficient,  Cq,  is  a  corrective  ftctor  to  bridge  the  gap  between 
the  ideal  and  real  worlds.  is  less  than  unity  due  to  the  presence  of 
the  orifice  boundary  layer  (viscosity  influence)  and  two  dimensional  flow 
effects . 
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As  the  gas  passes  through  the  orifice,  an  expansion  fan  is  formed  on 

the  inner  edge  of  the  orifice  and  spreads  out  into  the  flow  field.  The 

expansion  fan  is  initially  describable  in  terms  of  Prandtl-Meyer  theory, 
but  is  altered  from  that  form  as  it  proceeds  away  from  the  orifice  by 

axisymmetric  (as  opposed  to  planar)  flow  effects. 

The  initial  portion  of  the  supersonic  expansion  is  well  suited  for 
calculation  via  the  method  of  characteristics  (MOC),3  although  the 
interaction  of  the  boundary  layer  with  the  expansion  fan  is  difficult  to 
treat  analytically.  As  the  flow  proceeds  away  from  the  orifice  the 
rotational  and  randomized  translational  energy  is  transferred  to  directed 
kinetic  energy  of  the  flow.  The  flow  expansion  and  consequent  reduction  in 
collision  frequency,  causes  the  gas  to  depart  from  thermal  equilibrium.  A 
residual  rotational  energy  remains  in  the  gas  after  collisions  cease,  and 
this  is  often  characterized  by  a  rotational  temperature.  (The  populated 
rotational  states  do  not,  in  general,  adhere  to  a  Boltzmann  distribution, 4 
but  this  is  of  no  great  consequence  to  the  present  investigation.) 
Furthermore,  the  expansion  is  characterized  by  distinct  parallel  and 


!<■] 


3Shapiro,  A.H.,  The  Dynamics  and  Thermodynamics  of  Compressible  Fluid  Flow.  BW 

the  Ronald  Press  Co.,  New  York,  1953,  pp.  73-105.  v] 
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perpendicular  teaperatures  representing  residual  randoaized  kinetic  energy 
in  directions  parallel  to  and  perpendicular  to  the  Bean  flow  direction. 

The  portion  of  the  jet  at  large  angles  froa  the  symmetry  axis  is 
dominated  by  the  slower,  aore  easily  turned,  gas  froa  the  orifice  boundary 
layer.  Also,  the  lighter  aolecular  weight  species  become  aore  concentrated 
at  large  angles  froa  the  symmetry  axis  leaving  the  central  portion  of  the 
jet  aore  concentrated  in  species  with  heavier  aolecular  weights.5 

2.3  Real  Orifice  Flow  Effects 

In  principle,  it  would  be  possible  to  compute  the  flow  froa  the 
undisturbed  ataosphere  into  the  aass  spectrometer  in  a  single  solution.  In 
such  a  solution,  the  boundary  conditions  applied  would  correspond  to: 

1)  An  undisturbed  atmosphere  bounding  the  appropriate  portion  of 
the  solution  region. 

2)  A  wall  interaction  representing  the  molecular  collisions  with 
the  aass  spectrometer  casing  as  well  as  the  sides  of  the 
orifice. 

3)  Possible  downstream  conditions  if  the  solution  is  carried  to 
the  point  where  interaction  with  the  background  gas  within 
the  aass  spectrometer  is  important. 

Although  wall  boundary  conditions  do  have  a  degree  of  approxiaation 
that  reflects  our  incomplete  understanding  of  the  gas-surface  interaction 


4Sharafudinov,  R.G.  and  Skovorodko,  P.A.,  "Rotational  Level  Population 
Kinetics  in  Nitrogen  Free jets,"  Proceedings  of  the  12th  International 
Symposium  on  Rarefied  Gas  Dynamics,  AIAA  Press,  1980. 

5Bird,  G.A.,  "Breakdown  of  Continuum  Flow  in  Freejets  and  Rocket  Plumes." 
Proceedings  of  the  12th  International  Symposium  on  Rarefied  Gas  Dynamics, 
AIAA  Press,  1980. 
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process,  the  above  boundary  conditions  are  relatively  well  characterized. 
However,  the  above  procedure  does  involve  a  substantial  amount  of 
computational  effort  aimed  at  the  flow  outside  of  the  mass  spectrometer, 
since  that  flow  is  usually  characterized  by  a  relatively  small  mean  free 
path.  Since  the  primary  interest  of  the  present  study  is  to  describe  the 
flow  within  the  mass  spectrometer,  it  is  natural  to  try  to  reduce  the  size 
of  the  solution  region  while  retaining  the  essential  portion  of  the  flow. 
The  obvious  way  to  do  this  is  to  approximate  the  flow  at  the  orifice,  and 
merely  compute  the  flow  within  the  instrument. 

There  are  several  reasons  to  believe  that  such  a  procedure  would  be 
fruitful.  If  the  orifice  were  a  smooth  walled  converging-diverging  nozzle, 
then  the  flow  would  expand  to  a  Mach  number  of  unity  in  the  orifice  plane, 
and  it  would  be  known  exactly  (except  for  boundary  layer  and  two 
dimensional  effects)  irrespective  of  the  flow  downstream  of  the  orifice. 
Although  the  actual  orifice  has  a  sharp  edge,  intuition  implies  that  the 
downstream  flow  can  have  at  best  a  small  effect  on  the  orifice  flow,  so  it 
would  seem  plausible  to  specify  the  orifice  flow  a  priori.  Furthermore, 
the  present  study  is  concerned  mainly  with  the  flow  along  the  centerline  of 
the  resulting  jet,  where  boundary  layer  and  two  dimensional  effects  are 
minimized. 

Although  it  would  seem  to  be  a  classic  problem,  the  expansion  from  a 
uniform  state  through  a  sharp  edged  orifice  into  a  vacuum  has  apparently 
not  been  the  object  of  intense  study.  Liepmann6  considered  the  problem  in 
some  detail,  and  indicated  that  the  flow  in  the  plane  of  the  orifice  is 
subsonic  along  the  centerline,  going  supersonic  and  then  back  to  sonic  as 
the  edge  of  the  orifice  is  approached.  Little  else  was  said  in  this  or  any 

^Liepmann,  H.W.,  "Gaskinetics  and  Gasdynamics  of  Orifice  Flow,"  Journal  of 
Fluid  Mechanics,  Vol.  10,  No.  5,  1961,  pp.  65-79. 


8 


other  available  reference  about  the  actual  flow  in  the  orifice  plane. 
Measurements  were  presented,  but  they  were  limited  to  the  discharge 
coefficient,  which  was  subsequently  measured  by  Smetana,  Sherrill  and 
Schort7  to  greater  accuracy.  The  results  of  the  latter  study  are  presented 
in  Figure  3,  which  has  been  replotted  to  obtain  a  more  convenient  fora. 
(As  presented  in  Reference  7  both  axes  involved  the  unknown  mass  flow,  so 
an  iteration  was  required  to  determine  the  discharge  coefficient.  This 
problem  is  removed  when  the  Reynold's  number  is  defined  in  terms  of 
available  reference  quantities  rather  than  the  unknown  mass  flow.  The  two 
Reynold's  numbers  differ  only  by  a  factor  of  the  discharge  coefficient.) 

Use  of  the  discharge  coefficient  presented  in  Figure  3  accounts  for 
the  major  effect  of  a  sharp  edged  orifice  as  opposed  to  a  smooth  nozzle, 
and  it  is  apparently  the  only  available  effect  for  which  measurements 
exist.  The  approach  used  in  the  EXPANDO  code  is  to  adjust  the  orifice  mass 
flow  to  reflect  the  results  of  Figure  3,  but  the  mean  flow  velocity  is 
still  assumed  to  start  at  sonic  conditions.  (Individual  molecules,  of 
course,  will  have  a  distribution  of  velocities  about  the  Dean.)  Any  error 
incurred  via  this  procedure  should  have  little  effect  on  the  centerline 
solution  of  interest. 

2.4  Justification  of  Calculational  Approach 

For  the  most  part  the  departure  of  the  freejet  from  local  thermo¬ 
dynamic  equilibrium  marks  the  end  of  the  region  of  validity  for  continuum 
techniques  such  as  the  method  of  characteristics.  (An  exception  to  this  is 
a  recent  paper  by  Labowski,  Ryali,  and  Fenn  using  the  so-called 

7Smetana,  F.O.,  Sherrill,  W.A.  ,  and  Schort,  D.R.,  "Measurements  of  the 
Discharge  Characteristics  of  Sharp-edged  and  Round-edged  Orifices  in  the 
Transition  Regime,"  Proceedings  of  the  5'th  Rarefied  Gas  Dynamics 
Symposium,  Vol.  2.  Academic  Press,  1967. 
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Figure  3.  The  Discarge  Coefficient  for  a  Sharp  Edged  Orifice  as  Measured 
by  Smetana,  Sherrill  and  Schort.7  It  is  plotted  here  as  a  function  of 
Reynold's  number  based  on  sonic  conditions  and  the  orifice  diameter. 


nonequilibrium  method  of  characteristics  on  a  rotational ly  relaxing 
freejet.®  However,  the  method  does  not  apply  when  the  translational  modes 
go  out  of  equilibrium.  Since  the  region  between  rotational  and 
translational  nonequilibrium  is  usually  small,  the  method  is  inappropriate 
in  the  present  study.) 

The  only  presently  available  computational  technique  which  is  capable 
of  describing  the  critical  molecular  processes  occurring  in  the  latter 
portio: s  of  the  expansion  is  the  direct  simulation  Monte  Carlo  Method. 
There  was  no  choice  to  be  made  for  the  calculational  technique  to  be  used 
in  the  nonequilibrium  portion  of  the  expansion;  the  only  question  is 
whether  the  Monte  Carlo  method  should  have  been  used  for  the  entire 
solution  region,  or  whether  it  should  have  been  wedded  to  a  MOC  calculation 
of  the  equilibrium  portion.  The  answer  to  this  question  depended  largely 
on  how  large  the  equilibrium  region  was  expected  to  be.  (It  should  be 
stressed  that  the  Monte  Carlo  flow  field  is  valid  for  equilibrium  continuum 
flow  —  but  it  is  less  efficient  computationally.)  The  breakdown  of 
translational  and  rotational  equilibrium  has  been  considered  in  some  detail 
by  Bird,5,9  and  based  on  this  work  the  conclusion  is  that  a  MOC  calculation 
could  not  be  used  beyond  two  orifice  radii  from  the  orifice  at  an  altitude 
of  20  km  and  essentially  could  not  be  used  at  all  at  40  km  altitude.  This 
very  limited  range  of  validity  for  a  MOC  calculation  Indicated  that  the 
preferred  calculational  approach  was  simply  to  use  a  Monte  Carlo  method  for 
the  entire  flow  field. 


®Labowskl,  M.  Ryali.  S.,  and  Fenn,  J.B.,  "Flowfield  Calculations  in 
Nonequilibrium  Freejets  by  the  Method  of  Characteristics,"  Proceedings  of 
the  12th  International  Symposium  on  Rarefied  Gas  Dynamics,  AIAA  Press, 
1980. 

9Bird,  G.A.,  "Breakdown  of  Translational  and  Rotational  Equilibrium  in 
Gaseous  Expansions,"  AIAA  Journal.  Vol.  8,  No.  13,  November  1970.  pp. 
1997-2003. 
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3 .  CALCULATIONAL  GRID 


3.1  General  Considerations 

The  selection  of  grid  geometries  for  fluid  mechanic  calculations  can 
generally  be  regarded  as  more  of  an  art  than  a  science.  Considerations  in 
the  selection  of  a  grid  are: 

*  The  grid  should  be  as  simple  as  possible,  since  the  program 
must  repeatedly  decide  what  cell  molecules  are  in  as  they 
move.  If  this  required  the  solution  of  a  complex  equation, 
the  entire  program  would  run  significantly  slower  than  if  the 
cell  can  be  determined  easily. 

*  The  grid  should  concentrate  cells  where  gradients  are  the 

largest,  so  that  the  least  number  of  total  cells  (and 
molecules)  are  needed  to  obtain  an  accurate  solution. 

*  The  grid  should  provide  output  where  it  is  required,  with  the 
resolution  that  is  desired  for  the  answer  of  interest. 

The  EXPANDO  cell  structure  is  illustrated  schematically  in  Figure  4, 
and  is  characterized  by  the  following  relations: 

*  The  orifice  radius  is  divided  up  evenly  into  a  specified 

number  of  grid  spacings.  This  number  can  be  varied  to  achieve 
convergence . 

*  Above  the  orifice,  at  the  first  axial  location,  there  are  some 

additional  cells  of  the  same  size  to  help  describe  the 

expansion  that  occurs  around  the  orifice  edge.  The  number  of 
cells  above  the  orifice  can  also  be  varied  to  achieve 
convergence  along  the  flow  centerline  (the  region  of 
importance  for  the  present  investigation) . 

*  The  cells  are  constructed  between  a  series  of  planes  which  are 
perpendicular  to  the  jet  axis.  The  spacing  between  successive 
planes  increases  as  they  proceed  further  from  the  orifice, 
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Figure  4.  A  Drawing  of  a  Ceil  Mesh  of  the  Type  Used  in  Expando. 


since  the  flow  gradients  will  decrease  in  this  direction.  The 
rule  for  the  increase  of  the  spacing  between  the  planes  is 
given  by  the  equation: 


- hr-  ■  A  -  constant  .  (2) 

ro  +  zj 

where  Zj  is  the  location  of  the  J'th  plane  and  AZj  is  the 
spacing  between  it  and  the  previous  plane  (i.e.,  AZj  -  Zj  - 
Zj_j),  and  r0  is  the  orifice  radius.  This  relation  can  be 
solved  to  give  the  location  of  Zj  directly  as: 


Zj  '  rfl  [  (l-A)-J  -  1] 


(3) 


This  relation  is  easily  inverted  to  find  the  axial  cell 
location  of  a  particular  aolecule. 

The  cells  are  selected  to  have  a  radial  increment  equal  to  the 
axial  increment,  and  the  number  of  cells  per  axial  location  is 
kept  constant.  Hence,  as  the  cells  grow  larger  in  the  axial 
direction,  they  autoaatically  keep  pace  in  the  radial 
direction. 


3.2  Spatial  Segaentation  Of  The  Solution  Region 

The  code  was  written  so  that  it  has  the  ability  to  coapute  sequential 
spatial  segaents  of  the  solution  starting  froa  the  orifice.  (This  ability 
is  aerely  an  option  and  in  no  way  affects  the  capability  of  handling  the 
entire  flow  field  at  once.)  The  use  of  this  option  offers  the  potential  for 
a  substantial  decrease  in  the  aemory  and  computing  time  required  to  solve  a 
problem.  Details  of  the  spatial  segaentation  scheme  are  discussed  in  the 
subsections  below. 


3.2.1  Justification  for  Segaentation 

The  first  question  that  must  be  addressed  is  whether  the  segmentation 
of  the  solution  is  physically  and  mathematically  justified,  and  this 


question  is  critically  related  to  the  question  of  boundary  conditions.  The 
physical  laws  that  are  embodied  in  the  Monte  Carlo  solution  procedure, 
involving  molecular  translations  and  collisions,  are  as  valid  for  a  portion 
of  the  solution  region  as  they  are  for  the  whole  region.  In  order  to  carry 
out  the  solution  in  just  a  subregion,  however,  it  is  necessary  that  the 
boundary  conditions  can  be  specified  a  priori  along  the  boundaries  of  the 
subregion  in  question. 

For  a  Monte  Carlo  flow  field  calculation,  the  boundary  conditions  are 
imposed  by  specifying  the  velocity  distribution  function  for  incoming 
molecules  along  all  boundaries,  and  then  selecting  molecules  from  this 
distribution  with  the  proper  frequency  and  introducing  then  into  the 
simulation.  Usually  this  involves  extending  the  boundaries  to  a  region  of 
undisturbed  (known)  flow,  or  to  a  region  where  molecular  backflow  into  the 
solution  region  is  insignificant. 

In  the  simulation  of  mass  spectrometer  flow,  the  upstream  boundary  is 
taken  to  correspond  to  one  dimensional  sonic  conditions  at  the  orifice, 
except  that  the  mass  flow  is  reduced  by  an  empirically  determined  discharge 
coefficient.  The  solution  region  is  extended  far  enough  downstream  and  to 
the  side  so  that  the  flow  of  molecules  into  the  solution  region  from  these 
other  boundaries  can  be  neglected.  (An  exception  to  this  is  the  incursion 
of  background  gas  into  the  jet  which  could  be  treated  as  an  equilibrium  gas 
at  the  side  boundaries.  This  feature  is  not  included  in  the  code  nor  does 
it  affect  the  present  discussion  since  the  distribution  function  for  this 
gas  would  be  assumed  known  at  the  side  boundaries.)  Hence,  a  well  defined 
solution  can  be  carried  out  downstream  of  the  orifice  as  long  as  the 
downstream  boundary  is  far  enough  from  the  orifice  to  assure  that  molecular 
backflow  is  negligible.  Since  molecules  have  a  thermal  velocity  which  is 
on  the  order  of  the  speed  of  sound,  backflow  becomes  negligible  when  the 
local  Mach  number  is  large  compared  to  unity.  (It  is  interesting  to 
contrast  this  to  a  continuum  calculation,  which  requires  only  that  the 

1  5 


local  Mach  number  be  greater  than  one,  perhaps  by  a  very  small  amount,  lor 
there  to  be  no  upstream  influence.)  In  a  Maxwellian  gas,  the  probability 
of  an  individual  molecule  laving  an  upstream  velocity  direction  is  given  by 
P.  where 


P  =  -^erf  c  (M-vY/2  ) 


and  M  and  Y  denote  Mach  number  and  ratio  of  specific  heats,  respectively. 
For  sonic  flow  this  probability  is  on  the  order  of  5\,  but  by  the  time  the 
Mach  number  has  become  2.0  the  probability  is  on  the  order  of  4X10'-1  and, 
for  all  practical  purposes,  backflow  can  be  ignored. 

Hence,  the  imposition  of  the  downstream  boundary  condition  neglecting 
backflow  into  the  solution  region  is  justified  very  shortly  alter  the 
orifice,  since  the  flow  rapidly  becomes  substantially  supersonic.  As  long 
as  this  is  the  case,  it  is  perfectly  proper  to  solve  for  a  small  segment  of 
the  solution  region.  Once  that  solution  is  completed,  then  the  derived 
velocity  distribution  on  the  downstream  boundary  defines  the  required 
upstream  boundary  condition  for  the  next  solution  segment.  This 
information  is  automatically  written  to  a  file  which  is.  in  turn, 
automatically  read  as  input  for  the  next  segment. 


3.2.2  Advantages  of  Segmentation 


Although  the  physical  justification  for  segmenting  the  solution  region 
has  been  demonstrated,  there  remains  the  question  as  to  why  one  would  want 
to  "turn  one  problem  into  two  or  more  problems".  Frequently  the 
subdivision  of  a  physical  problem,  into  multiple  subproblems  reduces  the 
total  effort  involved,  and  this  is  no  exception.  Specific  advantages  at e 
< ■numerated  below. 
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3.2.2. 1  Reduced  Storage 


When  the  problem  Is  broken  up  into  segments,  the  total  nuaber  of  cells 
and  molecules  required  in  any  one  segment  is  less  than  required  for  the 
larger  single  solution.  This  translates  directly  into  a  decreased 
requirement  for  computer  core. 

3. 2. 2. 2  Decreased  Time  to  Achieve  Steady  Flow 

The  Monte  Carlo  technique  is  inherently  a  procedure  which  solves  an 
unsteady  flow  problem.  For  cases  such  as  the  present  one  where  the  problem 
of  interest  is  really  a  steady  state  flow,  this  is  solved  by  letting  the 
unsteady  solution  relax  to  a  steady  state.  The  computation  time  required 
to  achieve  a  steady  state  can  be  regarded  as  "computational  overhead"  for 
the  present  problem,  since  useful  steady  state  sampling  of  the  solution 
cannot  begin  until  steady  state  has  been  achieved. 

The  time  required  to  establish  a  steady  flow  can  be  estimated  a 
priori  as  the  length  of  the  solution  region  divided  by  the  orifice  flow 
velocity  (times  a  safety  factor).  Hence,  the  longer  the  initial  solution 
segment,  the  greater  is  the  period  of  unsteady  flow.  If  the  entire 
solvticr.  region  is  calculated  at  once,  then  the  program  may  well  spend  most 
of  its  effort  in  simply  achieving  steady  state.  If  the  solution  region  is 
segmented,  however,  then  the  first  segment  can  reach  steady  state 
substantially  faster  than  the  whole  solution  region  does.  This  is  of 
particular  significance  since  the  solution  near  the  orifice  is  the  most 
collision  dominated,  requiring  a  large  portion  of  the  computational  effort 
needed  in  each  time  step.  When  subsequent  segments  are  solved,  they  still 
require  a  relatively  long  time  to  achieve  steady  state  (though  the  time  is 
somewhat  diminished  since  the  solution  region  is  shorter),  but  the 
computational  effort  per  time  step  is  substantially  less 
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3. 2. 2. 3  Increased  Time  Step  for  Latter  Segments 

The  tine  step  in  a  Monte  Carlo  simulation  should  generally  be  small 
compared  to  the  mean  time  between  collisions  for  a  molecule,  since  the 
processes  of  translation  and  collisions  are  treated  separately  in  each  time 
step.  If  the  entire  solution  is  solved  for  at  once,  this  implies  that  the 
entire  solution  is  constrained  to  the  relatively  small  time  step  required 
by  the  collisional  region  near  the  orifice.  If  that  region  is  solved  for 
separately,  then  subsequent  segments  can  have  a  substantially  larger  time 
step  since  the  mean  time  between  collisions  is  larger  in  subsequent 
segments.  Hence,  even  though  the  latter  regions  still  require  a  relatively 
large  amount  of  simulation  time  to  achieve  steady  state,  this  time  is  more 
easily  accomplished  since  the  allowable  time  step  is  much  larger. 

3.2.3  Segmentation  Summary 

These  considerations  strongly  suggest  that  the  first  spatial  segment 
should  be  made  within  a  few  diameters  of  the  orifice,  since:  1)  The  flow 
is  by  then  already  sufficiently  supersonic  to  justify  the  neglect  of 
backflowing  molecules,  and  2)  The  number  of  time  steps  required  to  achieve 
steady  state  is  small,  which  is  particularly  important  in  this  collision 
dominated  region  of  the  flow.  Note  that  the  question  is  completely 
unrelated  to  whether  or  not  the  flow  is  in  equilibrium  (as  would  be  the 
case  if  an  initial  region  were  to  be  calculated  by  the  method  of 
characteristics).  It  is  simply  a  matter  of  making  the  calculation  of  the 
flow  field  more  efficient  by  separating  the  relatively  long  relaxation  time 
which  is  required  by  the  latter  portions  of  the  flow  field  from  the 
collision  dominance  which  is  characteristic  of  the  initial  portion.  Once 
the  first  segment  is  calculated,  then  subsequent  segments  benefit 
computationally  from  the  ability  to  take  substantially  larger  time  steps  in 
those  regions. 


GAS  MODEL 


4  . 


The  details  of  the  Monte  Carlo  solution  procedure  are  described 
extensively  in  References  1  and  2;  a  summary  of  the  majot  features  of  the 
method  is  provided  here  for  completeness.  Unlike  early  Monte  Carlo  models, 
the  molecules  are  allowed  to  possess  internal  energy,  and  interact  with 
velocity  dependent  cross  sections. 


4.1  Collision  Cross  Sections 


The  available  translational  eollisional  energy  between  the  two 
molecules,  Ec ,  is  given  in  terms  of  their  relative  velocity,  cr .  by 


E 


c 
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I^iicr 


(5) 


where  p1j  is  the  reduced  mass  of  the  pair;  i.e. 


1,1  Pi¬ 
rn  j  -n\j 


(h) 


and  mj  and  mj  represent  the  molecular  masses.  The  code  utilizes  the 
Variable-Hard-Sphere  (VHS)  model,  ®  where  molecules  have  a  collision  cross 
section  which  varies  as  an  inverse  power  of  the  available  collision  enejuy 
Hence,  if  ffjj  is  the  collision  cross  section  for  collisions  of  species  i 
with  species  j,  then  0A  j  is  given  by  a  relation  of  the  form 


^Bird.  G.  A.,  "Monte- Car lo  Simulation  in  an  Engineering  Context  . 
Proceedings  of  the  12th  International  Symposium  on  Rarefied  Gas  Dynan  us. 
Vol.  74,  Progress  in  Astronautics  and  Aeronautics.  A  ’  A  A ,  New  Vo:k  .  19tsl  . 


where  is  a  constant  coefficient. 

The  parameter  o  can  be  related  to  i),  the  exponent  of  distance  in  an 
inverse  power  intermolecular  force  law  via  the  relation1® 


Hence,  hard  sphere  molecules  (for  which  h  goes  to  infinity)  are  represented 
by  u  equal  to  zero.  There  is  a  substantial  body  of  evidence,  however,  that 
the  effective  size  of  molecules  does  indeed  decrease  with  increasing 
collision  energy,  so  a  positive  value  of  u  is  usually  a  better  choice,  a 
can  be  determined  from  molecular  beam  data,  or  from  its  macroscopic 
implications.  For  example,  if  s  is  the  exponent  for  the  variation  of  the 
viscosity  coefficient  with  temperature,  then  it  can  be  shown10  that 


=  u  +  0.5 


so  a  measurement  of  the  temperature  dependence  of  the  viscosity  coefficient 
serves  as  an  indirect  determination  of  o.  It  was  on  this  basis  that  a 
value  of  o  equal  to  approximately  0.24  is  recommended  for  the  major  species 
(air)  flow. 


Although  the  sizes  of  molecules  are  allowed  to  vary  in  the  VHS  model 
in  deciding  whether  or  not  a  collision  is  to  occur,  when  a  collision  does 
occur  the  post  collision  velocity  components  are  computed  as  if  it  were  a 
hard  sphere  collision.  This  results  in  a  substantial  computational 
simplification  and  yet  retains  good  agreement  with  the  macroscopic 


predictions  of  the  more  exact  model. 


4.2  Weighting  Factors 


Statistical  weighting  factors  are  a  crucial  element  of  a  successful 
Monte  Carlo  simulation,  allowing  trace  species  to  be  described  with 
reasonable  accuracy.  The  weighting  factor  is  the  number  of  "real" 
molecules  that  correspond  to  each  "simulated"  molecule.  A  "simulated" 
molecule  corresponds  to  one  molecule's  worth  of  storage  allocated  in  the 
program,  and  the  weighting  factor  is  its  statistical  weight.  The  weighting 
factors  used  in  EXPAXDO  are  dependent  on  cell  and  species.  The  weighting 
factors  are  dynamically  adjusted  as  the  simulation  proceeds,  keeping  the 
number  of  simulated  molecules  more  or  less  constant  while  allowing  the 
number  of  real  molecules  to  adjust  as  the  solution  evolves. 

4.3  Collision  Mechanics 

4.3.1  Elastic  Collisions 

Conservation  of  momentum  implies  that  the  center-of-a iss  velocity  of  a 
collision  pair  is  unchanged  by  the  collision;  and  conservation  of  energy 
then  implies  that  the  magnitude  of  the  relative  velocity  between  the 
collision  partners  is  also  unchanged  by  the  collision.11  Since  the 
collision  is  treated  as  a  statistical  event,  all  that  remains  is  to  select 
the  direction  of  the  post  collision  relative  velocity  vector  from  the 
correct  distribution.  As  mentioned  earlier,  collisions  in  the  VHS  mouel 
are  treated  as  have  sphere  collisions  when  they  occur  (though  they  cio  not 
occur  with  the  same  velocity  dependence  as  do  hard  sphere  collisions).  For 
such  molecules,  ail  directions  for  the  post  collision  relative  velocity 

^Viment.,  Walter  G.,  and  Kruger,  Charles  H.,  Jr.,  1  nt  rodui  t  ;  cm  '  <• 
Physjca,  Gas  Dynamics,  John  Wiley  and  Sons.  19(>r> .  pp .  3  5C-35C 
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vector  are  equally  likely.  This  is  the  chief  computational  simplicity  of 
the  VHS  model.  An  elastic  collision  is  simulated  by  selecting  a  new 
direction  for  the  relative  velocity  vector  and  then  computing  the 
corresponding  molecular  velocities. 

4.3.2  Inelastic  Collisions 

Molecules  are  allowed  to  have  internal  energy  which  is  described  via 
the  phenomenological  model  of  Borgnakke  and  Larsen.12  jn  this  model, 
transfer  of  energy  between  internal  and  translational  modes  is  allowed,  but 
it  is  necessary  to  assume  that  each  species  has  a  fixed  number  of  internal 
degrees  of  freedom,  Cj.  A  collision  is  assumed  to  be  either  perfectly 
elastic  or  perfectly  inelastic  via  a  user  specified  probability.  A 
perfectly  inelastic  collision  is  achieved  by  summing  the  total 
pre-collision  energy  (i.e.,  the  internal  energy  of  both  molecules  plus  the 
translational  energy  of  their  relative  motion)  and  then  assigning  post 
collision  values  from  the  equilibrium  distribution  for  collisions  with  that 
total  amount  of  energy,  taking  into  account  the  number  of  internal  degrees 
of  freedom  in  the  two  molecules.  Note  that  this  model  has  the  ability  to 
relax  from  a  nonequl librium  to  an  equilibrium  state  via  an  effective 
collision  number.  The  ability  to  exchange  internal  energy  in  such  a  manner 
comprises  a  significant  increase  in  capability  for  Monte  Carlo  codes  beyond 
the  previous  models  where  molecules  had  no  internal  energy.  It  is  this 
capability  which  enables  the  codes  to  realistically  predict  the  macroscopic 
effects  of  polyatomic  gas-  flow.  This  model  predicts  an  overall  level  of 
internal  (usually  rotational)  energy,  but  it  does  not  describe  individual 


12Borgnakke,  Claus,  and  Larsen,  Paul  S.,  "Statistical  Collision  Model  for 
Monte  Carlo  Simulation  of  Polyatomic  Gas  Mixture",  Journal  of  Computational 
Physics,  Vol .  18.  1975,  pp.  405-420. 
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quantum  states.  When  this  is  desired,  it  can  be  accomplished  by 
considering  the  separate  states  to  be  distinct  species. 

Let  Cj  and  C2  be  the  nuaber  of  internal  degrees  of  freedoa  of  the  two 
molecules  in  an  Inelastic  collision,  and  Eg  be  the  total  collision  energy 
defined  by 

Es  *  Eci  +  Eli  +  E2i  •  d°) 

where  Ecj  is  the  initial  translational  collision  energy  defined  by  Eq.  (5), 
and  Ej^  and  E2i  are  the  pre-collision  internal  energies  of  the  two 
molecules.  The  somewhat  cumbersome  expressions  given  in  Ref.  12,  can  be 
recast  in  terms  of  the  Chi-Square  distribution.  Post-collision  values  for 
the  respective  energies  are  given  by 

xlEs 

Eif  XX  *  X2  +  X3  1  (11) 


x2es 

E2f  *  Xj  +  X2  +  X3 

and 


(12) 


X3Es 

Ecf  =  Xj  +  X2  +  X3  ’  (13) 

where  Xj  is  selected  from  a  Chi-Square  distribution  with  Sj  degrees  of 
freedom,  X2  is  selected  from  a  Chi-Square  distribution  with  S2  degrees  of 
freedom  and  X3  is  selected  from  a  Chi-Square  distribution  with  2(2~o) 
degrees  of  freedom.  The  post  collision  translational  energy  is  then  used 
to  determine  a  new  relative  velocity  between  the  two  molecules. 
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4.3.3  Reactive  Collisions 

Reactive  collisions  can  be  simulated  directly  in  the  EXPANDO  code. 
The  treatment  of  reactive  collisions  is  similar  to  that  for  inelastic 
collisions,  except  that  a  heat  of  reaction  is  added  to  the  total  energy 
expressed  in  Eq.  (10).  The  reaction  is  simulated  with  a  probability  which 
is  proportional  to  the  ratio  of  the  reactive  to  collision  cross  section  at 
the  relative  velocity  for  the  collision.  Reactive  collisions  can  result  in 
the  disappearance  of  reactant  molecules,  with  the  post  collision  state 
being  applied  to  the  product  molecules. 

5.  KINETICS  OP  PROTON  HYDRATE  REACTIONS 

Laboratory  and  in  situ  mass  spectrometer  measurements13  indicate  that 
stratospheric  positive  ions  undergo  significant  dissociation  in  the 
presence  of  electric  fields  of  moderate  strengths  (0  ■*  20  V/cm).  In  order 
to  understand  this  type  of  phenomena,  we  have  investigated  the  kinetics  of 
the  proton  hydrate  (PH)  reactions.  The  PH  system  was  chosen  for  two 
reasons:  (1)  in  situ  measurements  of  low  hydration  order  PH  show 
abundances  much  greater  than  theoretical  (thermodynamic)  predictions  and 
(2)  there  are  numerous  experimental  results  on  the  kinetics  and 
thermodynamics  of  PH  reactions.  In  the  remainder  of  this  section,  we 
discuss  the  experimental  results  pertinent  to  the  PH  reactions,  as  well  as 
the  approximate  theories  for  the  agglomeration/fragmentation  processes, 
which  are  used  in  the  Monte  Carlo  simulation. 


13Arnold,  F.  ,  Henschen,  G.  ,  and  Ferguson,  E.  E.  ,  "Mass  Spectrometric 
Measurements  of  Fractional  Ion  Abundances  in  the  Stratosphere-Positive 
Ions,”  Planet.  Space  Sci.  29(3983),  pp.  385-193. 
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5.1  Laboratory  Measureaents  of  PH  Reactions 

The  overall  reaction  in  which  we  are  interested  is 

kf 

H20+H+(H20)n+M  *  H^(H20)n+14M  ,  (14) 

kr 

where  n=l,2 . 5  and  M  represents  the  third  body.  Usually,  such  ion 

agglomeration/fragmentation  reactions  are  represented  by  the  following 
mechanism: 


kc 

H20+H+(H20)n  -»  [H+(H20)n+1]*  .  (15) 

ks 

M+[H+  (H20)n+1]*  -»  M+H+  (H20)n+1  .  (16) 


*  kd 

[H+(H20)n+1]  H20+H+(H20)n 

and 

ka  , 

M+H+(H20)n+1  -*  M+[H+(H20)n+ir 


(17) 


(18) 


[For  notational  simplicity,  the  cluster  size  index,  n,  has  been  suppressed 
in  writing  the  rate  constants.  It  should  be  apparent  to  the  reader  that 
kf ,  kr,  etc.  refer  to  the  appropriate  rate  constant  for  a  particular 
cluster  size,  i.e.,  kfn,  hr  n+1,  etc.].  Assuming  steady  state  for 
[H+(H20)n+j]*  allows  the  overall  forward  and  reverse  rates  to  be  written  in 
terms  of  the  Individual  rates  for  Eqs.  (15)-(18)  as 


kcksW 
f  =  kd+ks[M] 


(19) 


♦ 
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In  the  low  pressure  Halt  (l.e..  k(j>>ks[M])l  Eqs.  (19)  and  (20)  reduce  to 


kf 


kcks 


(21) 


and 


kr  -  ka[M]  -  kJ[M]  (22) 

There  have  been  extensive  experiaental  studies  to  determine  the  rate 

constants,  k°,  as  a  function  of  PH  cluster  size  with  M=N2,14  0215  an<^ 

CH4.16  For  M=CH4,  the  temperature  dependence  of  k^  was  also  determined  over 

a  wide  temperature  range  for  n=l,2 . 5.  The  latest  (and  probably  most 

accurate)  measurements  of  the  equilibrium  constant  (K=kf/kr)  are  those  of 

Kebarle  et  al.16  It  is  important  to  note  that  there  are  no  direct 

measurements  of  the  reverse  rate,  k°  However,  values  of  k°  have  been 

r  r 


14Good,  A.,  Durden,  D.  A.,  and  Kebarle,  P.,  "Ion  Molecule  Reactions  in  Pure 
Nitrogen  and  Nitrogen  Containing  Traces  of  Water  at  Total  Pressures  0.5-4 
torr.  Kinetics  of  Clustering  Reactions  Forming  H*(H20)n,"  J.  Chem.  Phys. 
52(1970),  pp.  212-221. 

15Good,  A.,  Durden,  D.  A.,  and  Kebarle,  P.,  "Mechanism  and  Rate  Constants 
of  Ion-Molecule  Reactions  Leading  to  Formation  of  H+(H20)n  in  Moist  Oxygen 
and  Air,”  J.  Chem.  Phys.  52(1970),  pp.  222-229. 

16Lau,  Y.  K.,  Ikuta,  S.,  and  Kebarle,  P.,  "Thermodynamics  and  Kinetics  of 
the  Gas-Phase  Reactions:  Ho0+(H?0)n_1+Ho0?*Ho0+(Ho0)n, "  J.  Am.  Chem.  Soc. 
104(1982).  pp.  1463-1469. 


obtained  using  the  forward  rates  and  equilibrium  constants.  The 
experimental  data  which  are  most  relevant  to  the  present  work  are 
summarized  in  Table  1. 


Table  1.  Properties  of  the  N2+H20+H* (H20)n  *  N2+H+ (H20)n+1 

Reactions. 


n  Kn,n+18 

(cm3) 

kf  ,n 

(cm®/sec) 

^r ,n+l 
(cm3/sec) 

-“Wia 

(kcal/mole) 

bd 

1  2.060X10"2 

3.4xl0-27b 

1 . 7xl0-25 

31.6 

4 

2  1.207xl0"10 

2.3xl0-27b 

1 .9xl0~17 

19.5 

7.5 

3  3.405X10"10 

2.4xl0"27b 

7.0xl0'15 

17.9 

8.1 

4  4.909X10'16 

1.5x10*27c 

3.1xl0~12 

12.7 

14 

5  3.966xl0"17 

— - 

— 

11.6 

15.3 

“From  ref.  16 

bFrom  ref.  14 

cFrom  ref.  15 

Adjusted  to 

300  K  and  M*N2. 

dFrom  ref.  16 

The  parameter  b  gives  the 

temperature 

dependence  of  kf  i.e.,  kfn*T-t>. 

The  data  given  in  Table  1  will  provide  the  necessary  checks  and 
(sometime?  >  inputs  for  the  theoretical  cross  sections  which  will  be 
developed  in  the  following  sections. 


5.2  Rate  of  (H*(H20)n+1]*  Decomposition 

The  rate  of  unimolecular  decomposition  (Eq.  (17))  can  be  estimated 
from  classical  RRK  theory.1'  It  has  recently  been  shown1®  that  a  classical 
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RRK  description  of  Ar  clustering  provides  a  semi-quantitative  description 
of  the  uniaolecular  decomposition.  The  form  of  the  (1st  order)  rate 
constant  is 

kd  *  A((E-E0)/E)S"1  ,  (23) 

where  the  constant  A  is  the  limiting  rate,  E  is  the  total  energy  of  the 
complex,  E0  is  the  critical  energy  for  dissociation,  and  s  is  the  effective 
number  of  oscillators  in  the  complex.  In  the  present  case,  A  is  chosen 
arbitrarily  (albeit  judiciously),  E0  is  the  thermodynamic  dissociation 
energy,  and  s  is  approximately  equal  to  b+1,  where  b  is  the  temperature 
dependence  of  the  measured  forward  rate  (see  Table  1).  The  total  energy  of 
the  complex  is  given  by 

E  =  •|mv2+E0+sRT  ,  (24) 

where  we  have  assumed  the  internal  energy  of  the  complex  is  sRT.  Table  2 
contains  the  parameters  appropriate  to  the  PH  complexes.  It  is  interesting 
to  note  that  the  values  of  s  deduced  from  experiment  correlates  very  nicely 
with  3N-6  ( 3N-5  for  a  linear  complex),  if  one  assumes  that  the  H20  ligands 
have  no  internal  degrees  of  freedom. 


*7See,  for  example,  Forst,  W.  ,  Theory  of  Unimolecular  Reactions  (Academic 
Press,  N.  Y. ,  1973) .  “  . 

18Brady,  J.  W. ,  Doll,  J.  D. ,  and  Thompson,  D.  L.,  J.  Chem .  Phys .  73(1980) , 
pp.  2767-2772. 
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Table  2.  Parameters  for  RRK  Treatment  of  Cluster  Dissociation. 


1 

2xl013 

2.1955xl0-12 

5 

4 

2 

lxlO14 

1.3548X10-12 

8 

6 

3 

lxlO14 

1.2437X10-12 

9 

9 

4 

lxlO14 

8.8237X10-13 

15 

12 

5 

lxlO14 

8.0595X10-13 

16 

15 

5.3  Theoretical  Cross  Sections  for  PH  Reactions 

The  cross  sections  for  the  above  processes  can  usua3'y  be  obtained  in 
analytical  form  using  simplified  theories  of  ion-molecule  collisions. 19,20 
In  general,  rate  constants  for  ion-molecule  reactions  display  either  no 
temperature  dependence  or  a  negative  temperature  dependence.  Since 

19Castleman,  A.  W.,  "Nucleation  and  Molecular  Clustering  About  Ions." 
Advances  in  Colloid  and  Interface  Science  K)(1979),  pp.  73-l?8.  and 
references  therein. 

20Hsieh,  E.  T.  -Y.  and  Castleman,  A.  W.,  "A  Reconsideration  of  the  Theory 
of  Capture  Cross  Sections  for  Ion/Molecule  Reactions  and  a  Total  Energy  and 
Angular  Momentum  Conserved  Average  Charge-Dipole  Interaction  Theory," 
International  Journal  of  Mass  Spectrometry  and  Ion  Physics  40(1981),  pp. 
295-329 . 


29 


I  .’b  X 


ion-aolecule  cross  sections  are  usually  controlled  by  long  range  forces, 
collisions  not  only  lack  an  activation  energy,  but  they  also  possess  large 
cross  sections.  The  first  seai-quantitative  treatment  of  ion-aolecule 
collisions  involved  the  use  of  the  Langevin  expression  for  an  ion-induced 
dipole  long  range  interaction.  Interestingly,  the  rate  constant  calculated 
using  the  Langevin  approach  shows  no  temperature  dependence.  Since  the 
Langevin  expression  shows  the  essential  features  of  the  collision  theories 
appropriate  for  the  present  study,  we  will  consider  this  approach  in 
detail.  The  extension  of  the  Langevin-type  treatment  of  collision  dynamics 
should  then  be  obvious. 


The  basic  idea  in  the  Langevin  approach  is  that  the  collision  cross 
section  is  completely  determined  by  long  range  interactions,  allowing  an 
extremely  complicated  N-body  scattering  problem  to  be  reduced  to  a  simple 
two-body  problem.  The  total  energy  for  the  two-body  problem  is  given  by 

{mv2  =  TR+Veff  ,  (25) 


where  p  is  the  reduced  mass  between  the  collision  partners,  v  is  the 
initial  relative  velocity,  TR  is  the  relative  nslational  energy.  Veff 
is  the  effective  potential  energy  given  by 


V 


eff 


(Mvb)2  ae2 
2pR2  2R4 


(26) 


where  b  is  the  impact  parameter,  a  is  the  polarizability  of  the  neutral 
molecule,  and  e  is  the  electronic  charge.  The  effective  potential,  Vefj. 
defined  by  Eq.  (26),  possesses  a  maximum  value  for  a  particular  value  of 
R=RC,  which  can  be  obtained  from 


dVeff(v.b,R) 


~4  -  ~W  -»  ■.->  ", 


For  collisions  with  a  relative  translational  energy,  TR.  satisfying 

Tr  >  MV2-Veff(v,b,Rc)  ,  (28) 

the  two  particles  have  enough  energy  to  surmount  the  centrifugal  barrier 
(»Veff (v,b,Rc))  and  enter  the  strong  collision  region.  If  Eq.  (28)  is  not 
satisfied,  then  the  particles  will  be  reflected  and  no  chemical  interaction 
can  occur.  In  Langevin  theory,  it  is  assumed  that  once  the  particles  enter 
the  strong  collision  region,  the  reaction  occurs  with  unit  probability. 
The  centrifugal  barrier  height  depends  on  both  v  and  b.  For  a  specified  v, 
the  barrier  height  will  increase  with  b,  until  an  impact  parameter  is 
reached  such  that 


|ivz  =  Veff(v,bc.Rc)  .  (29) 

Notice  that  in  Eq.  (29),  the  value  of  TR  is  taken  to  be  zero.  From  the 
above  discussion,  the  only  collisions  leading  to  reaction  are  those  with 
b<bc.  Therefore,  the  cross  section  for  ion-molecule  reactions  is  given  by 


The  value  of  bc  can  be  obtained  from  Eqs.  (26),  (27),  and  (29)  as 


2.1 


>c  =  — 'V4ae‘:/p 


Notice  that  the  expression  for  the  Langevin  cross  section  possesses  an 
inverse  velocity  dependence,  which  results  in  a  rate  constant  independent 
of  temperature.  This  behavior  is  easily  subsumed  within  the  Monte  Carlo 
formalism  by  choosing  y  =  0.5  in  Eq .  (7). 
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5.4  Complex  Activation  and  Stabilization  in  the  Monte  Carlo  Model 


The  ion-neutral  collision  cross  section,  a^,  is  obtained  from  Langevin 
theory.  There  is  no  ambient  temperature  dependence,  only  a  dependence  on 
the  polarizability  of  the  neutral  molecule  and  the  reduced  mass  of  the 
collision  partners.  The  collision  rates  are  given  in  Table  3  for  M=N2- 


Table  3.  Collision  Rate,  vffL,  for  N2+[H+(H20)n+1] . 


n 

vffL 

(cm3/sec) 

1 

7 . 7158x10”*° 

2 

7.1514X10"10 

3 

6.8474X10'10 

4 

6.6570X10-10 

5 

6.5265X10-10 

Within  the  Monte  Carlo  framework,  cluster  ions  are  treated  as  "big 
molecules”  with  2s  internal  degrees  of  freedom.  There  is  no  inherent 
difference  in  the  treatment  of  activated  and  stable  clusters  other  than  the 
possession  of  Internal  energy  greater  than  the  dissociation  energy.  When  a 
cluster  is  activated,  it  is  assigned  a  lifetime  from  an  exponential 
distribution  with  the  mean  value  determined  by  Eq.  (23).  Every  time  a 
cluster  (activated  or  not)  undergoes  a  collision,  it  is  treated  as  an 
inelastic  collision  (see  Section  4.3.2),  and  a  new  internal  energy  and 
lifetime  result.  Hence,  a  cluster  can  be  activated  or  stabilized  in  a 
collision  depending  upon  the  post-collision  internal  energy  of  the  cluster 
ion.  When  it's  lifetime  comes  up,  a  water  molecule  is  dissociated  off  and 
the  dissociation  energy  is  removed  from  the  product  pair. 


5.5  Cross  Section  for  Coaplex  Formation 


The  cross  section  for  coaplex  formation  (see  Eq.  (15)),  <jc,  cun  be 
obtained  from  the  average  dipole  orientation  (ADO)  theory.21  The  ADO  theory 
is  an  extension  of  Langevin  theory  which  allows  the  neutral  molecule  (i.e., 
H20)  to  possess  a  permanent,  as  well  as  induced,  dipole  aoaent.  In  this 
case,  the  effective  potential  is  given  by 


vef  f  = 


(Mvb)2 

2pR2 


Hde 

— — <cose(R)> 
R2 


(32) 


where  pd  is  the  dipole  moment  and  e(R)  is  the  dipole  orientation  angle. 
The  last  term  in  Eq.  (32)  represents  the  orientation  averaged  charge-dipole 
interaction.  Following  the  saae  procedure  outlined  above,  we  obtain  for  Rc 
and  bc 

2  oe2  ^de  d<cose>  , 

v  ‘  ^7  *  5C-B-  1331 

and 


b 


2 

c 


_2  ae*- 

R‘  +  — - — - 
c  pv2Rc2 


2Mde 

- — <COS0> 

pv2 


(34) 


Eqs.  (33)  and  (34)  are  implicit  relations  for  Rc  and  bc,  which  aust  be 
solved  numerically.  The  cross  sections  are  then  obtained  via  Eq.  (30). 


The  velocity  dependent  capture  rate  constants  vcc  for  H20+(H*H20)n 

(n=l,2 . 5)  are  shown  in  Fig.  5,  assuming  an  ambient  temperature  of  22? 

K.  The  (slight)  ambient  temperature  dependence  is  due  to  the  dipole 
rotation.  At  velocities  near  room  temperature,  the  rate  constant  is  nearly 


21Su,  T.  and  Bowers,  M.  T.,  J.  Chem.  Phys .  58^(1973),  pp.  3027-3037.  See 
Ref.  20  for  criticism  and  corrections  concerning  this  approach. 
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Figure  5.  ADO  Rate  Constants  for  Cluster  Formation  as  A  Function  of 
Cluster  Size  and  Collision  Velocity. 


temperature  independent.  As  the  collision  velocity  increases,  the  capture 
rate  decreases  as  1/v  until  the  Langevin  (high  velocity)  liait  is  reached. 
Notice  that  ADO  theory  predicts  a  rate  constant  up  to  an  order  of  magnitude 
greater  than  the  corresponding  Langevin  rate.  the  capture  rate  constant 
decreases  as  the  cluster  size  increases,  due  to  the  increase  in  the  reduced 
■ass  of  the  collision  partners. 


As  aentioned  above,  the  result  is  not  in  closed  form.  Therefore,  the 
calculated  cross  sections  are  fit  to  the  following  expression: 


vflrc  =  ^ 


(v<vn*) 


(35) 


and 


vor„  = 


bn 

T+Cn 


(y>vn*) 


(36) 


The  parameters  resulting  from  the  fit  of  Eqs.  (35)  and  (36)  to  the  rates 
shown  in  Fig.  5  are  given  in  Table  4. 


Table  4.  Parameters  Giving  the  Velocity  Dependence  of  the  Rate 
for  the  Process  H20+H+(H20)n  -»  lH+(H20)n+1]*. 


n  "n  un  ''n 

(ca/secj  (cm3/sec)  (cm4/sec2)  (cm3/sec) 


1 

4.9xl04 

8.316X10-9 

3.616X10'4 

9.372X10"10 

2 

4 . 3xl04 

7 . 261x1 0-9 

2.770X10"4 

8.188X10'10 

3 

4  .  lxlO4 

6.791xl0"9 

2 , 468xl0-4 

7 .721x10'*° 

4 

4 .  OxlO4 

6.572X10-9 

2 .329xl0*4 

7 . 500x10" 

5 

3.9xl04 

6.476X10-9 

2.239X10'4 

7 . 350xl0~10 

6.  SAMPLE  CASES  WITH  DISCUSSION 


The  EXPANDO  code  was  exercised  for  two  cases  to  produce  quantitative 
results.  These  cases  correspond  to  a  single  instrument  operated  at 
altitudes  of  30  and  38  kilometers,  and  were  designed  to  represent  actual 
experimental  cases  as  closely  as  possible.  These  cases  are  described  in 
detail  in  the  subsections  below. 


6.1  Instrument  Parameters 

The  instrument  was  characterized  in  the  code  by  the  parameters  given 
in  Table  5.  These  are  the  only  parameters  of  the  instrument  which  are  used 
by  the  EXPANDO  code,  and  they  are  felt  to  be  the  Important  ones  in 
determining  the  character  of  the  initial,  collision  dominated,  portion  of 
the  expansion  which  is  of  primary  concern  to  the  code. 


Table  5.  Instrument  Parameters  Input  to  EXPANDO  for  the  Sample 
Calculations. 


Parameter 

Value 

Units 

Orifice 

Diameter 

0.0406 

cm 

Orifice 

to  Skimmer  Distance 

2.54 

cm 

Skimmer 

Diameter 

0.9728 

cm 

Orifice 

Potential 

-1.0 

volt 

Skimmer 

Potential 

-50.0 

volt 

6.2  Atmospheric  Parameters 

The  1976  standard  atmosphere  was  used  to  obtain  the  pressure  and 
temperature  at  the  two  altitudes,  and  the  water  concentration  was  obtained 
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fro*  Ref.  22.  These  inputs  were  then  used  in  a  Model  given  in  Ref.  23  to 
predict  initial  ion  cluster  concentrations  to  use  in  starting  the 
calculation  at  the  orifice.24  These  parameters  are  sumarized  in  Table  6. 


Table  6.  Atmospheric  Parameters  Used  in  the  Sample  Calculations. 


Parameter 


Ambient  Number  Density  (cm-3) 
Ambient  Temperature  (K) 
H+(H2°)3  Mole  Fraction 
H+(H20)4  Mole  Fraction 
H+(H20)5  Mole  Fraction 
H*(H20)6  Mole  Fraction 
N2  Mole  Fraction 
02  Mole  Fraction 
H20  Mole  Fraction 


30  Km  38  Km 

Value  Value 


3.830xl017 

226.5 

3.43xl0~18 

1.81xl0~15 

3.27xl0~15 

1.41xl0~16 

0.7885 

0.2115 

4.0x10"® 


1. 116X1017 
244.8 
1.19X10"17 
6.20X10-15 
1.12X10"14 
4.84X10"16 
0.7885 
0.2115 
5.5x10"® 


The  velocity  dependence  of  the  cross  sections  for  N2,  02  and  H20  were 
determined  from  gas  viscosity  data,25  using  relations  given  in  Ref.  10,  and 
the  number  of  internal  degrees  of  freedom  was  determined  from  their  ratio 


22Deguchi,  S.  and  Muhleman,  D.  0.,  "Mesospheric  Water  Vapor" .Journal  of 
Geophysical  Research,  87(C2),  Feb.  20,  1982,  pp.  1343-1346. 

23Kebarle,  P.,  Annual  Review  of  Physical  Chemistry.  28,  1977,  p.  445. 

24Bc llanthin,  J.,  Private  Communication. 

25Handbook  of  Chemistry  and  Physics,  55th  Edition,  CRC  Press,  pp.  F58. 
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of  specific  heats.  (The  ratio  of  specific  heats  was  taken  at  150  K,  where 
the  rotational  Bodes  of  air  are  not  fully  excited.  This  led  to  less  than 
two  internal  degrees  of  freedoa  for  N2  an(*  ^2  ■ )  A  weighted  wean  value  of 
t»=0.2428  was  used  in  the  aajor  species  solution,  and  the  reference  values 
in  Table  7  then  define  the  variable  Ajj  in  Eq.  (7).  The  probability  of  a 
collision  being  elastic  was  taken  to  be  50%,  except  for  collisions 
involving  cluster  ions  which  were  always  taken  to  be  inelastic.  A  coaplete 
listing  of  the  cluster  ion  parameters  which  were  used  was  presented  in 
Section  5. 


Table  7.  Air  Molecular  Paraneters  Used  in  the  Saaple 
Calculations . 


Molecule 

Reference 
Cross  Section 
(cm2) 

Reference 

Velocity 

(ca/s) 

Energy 
Exponent ,u 

Internal 
Degrees  of 
Freedoa, C 

n2 

5.902xl0-15 

5.619X104 

0.2277 

1.444 

°2 

6.101X10-15 

5.149X104 

0.2992 

1.640 

h2o 

1.628xl0~14 

6.182xl04 

0.6211 

3.000 

6.3  Calculational  Results 

The  probability  of  a  aolecule  undergoing  a  reaction  per  unit  distance 
traveled  is  given  by  P,  where 

P  =  nr/(nvz)  ,  (37) 

and  nr  is  the  volumetric  reaction  rate  (reactions  per  unit  volume  per  unit 
tiae),  n  is  the  local  number  density  of  reactant  (molecules  per  unit 
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voluae),  and  v2  is  the  local  axial  velocity.  P  gives  a  direct  aeasure  of 
the  laportance  of  a  process  as  a  function  of  position  within  the  Bass 
spectrometer.  This  parameter  is  shown  in  Figure  6  for  the  fragmentation 
(uniaolecular  decomposition)  and  agglomeration  processes  acting  on  H'*' ( H20 ) 5 
at  30  km  altitude.  The  minor  wigiles  in  the  curves  are  a  result  of 
statistical  scatter  in  the  Monte  Carlo  solution  and  have  no  physical 
significance.  This  scatter  does  not  interfere  with  the  interpretation  of 
the  curves. 

The  fragmentation  process  does  not  start  until  approximately  0.2 
centimeters  from  the  orifice,  but  it  then  rises  rapidly  and  dies  off 
relatively  slowly  downstream.  This  behavior  is  readily  understood  in  terms 
of  the  relevant  physics.  Firstly,  the  uniaolecular  decomposition  of 
clusters  is  primarily  a  result  of  collisional  excitation  of  the  clusters. 
Hence,  it  depends  mainly  on  collision  frequency  and  collision  energy. 
These  quantities  are  both  considerable  different  for  ion-neutral  collisions 
than  they  are  for  neutral-neutral  collisions.  As  the  flow  expands,  the 
number  densities  decrease  substantial ly .  The  relative  velocity  between 
neutrals  drops  considerably  due  to  the  decreasing  temperature,  but  the 
relative  velocity  between  ions  and  neutrals  increases  due  to  acceleration 
of  the  ions  by  the  electric  field.  The  result  is  that  the  collision 
frequency  does  not  decrease  nearly  as  rapidly  for  ion-neutral  collisions  as 
for  neutral-neutral  collisions.  Similarly,  when  collisions  occur,  they 
occur  at  higher  energies.  All  these  factors  combine  to  produce  the 
relatively  gentle  slope  in  the  fragmentation  rate  curve. 

The  agglomeration  process  is  most  important  at  the  orifice,  and 
decreases  downstream  to  very  small  levels.  Although  agglomeration  is  also 
an  ion-neutral  process,  it  has  a  cross  section  which  decreases  with 
i  creasing  energy  (see  Fig.  5),  so  it  has  a  somewhat  greater  slope  in  the 
downstream  region.  However,  the  primary  issue  for  agglomeration  is  that  it 
involves  the  interaction  of  a  cluster  with  a  minor  species  (H20),  which  has 
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a  concentration  of  only  a  few  parts  per  Million.  There  are  not  enough 
collisions  between  the  clusters  and  water  for  agglomeration  to  be  a 
significant  process  in  these  calculations.  This  conclusion  is  made 
quantitatively  in  Pigs.  7  and  8  for  the  30  and  38  km  cases.  In  these 
figures,  the  integrated  probability  of  a  cluster  acquiring  another  water 
■olecule  is  presented  as  a  function  of  axial  distance.  At  30  ka  the 
probability  is  always  less  than  1%,  and  at  38  ka  its  aaxiaua  is  about  l^*. 

Figures  9  and  10  show  the  integrated  loss  from  fragmentation  for  the 
two  altitudes.  In  contrast  to  the  case  for  aggloaeration ,  it  can  see  that 
the  probability  of  a  cluster  being  fragmented  in  the  aass  spectroaeter  is 
on  the  order  of  50%  for  30  ka  and  25%  for  38  ka.  The  value  is  less  for  the 
higher  altitude  case  due  to  the  fewer  collisions  that  occur  at  lower 
density,  but  the  perturbation  of  the  Measured  values  would  have  to  be 
considered  significant  at  either  altitude. 

Figures  7-10  showed  the  loss  rate  for  a  particular  type  of  cluster, 
without  considering  gain  from  any  source.  In  a  Measurement,  of  course,  the 
interpretation  of  the  results  is  complicated  by  the  fact  that  when  a 
cluster  fragments  it  produces  another  cluster  with  one  less  water  molecule. 
The  calculated  ion  cluster  number  densities  are  plotted  as  a  function  of 
axial  position  in  Figures  11  and  12  for  the  two  altitudes.  The  most 
obvious  feature  of  these  calculations  is  the  abrupt  increase  in  the  number 
density  of  H+(H20)3  clusters  in  the  region  between  0.2  and  0.4  centimeters 
from  the  orifice.  This  feature  is  a  result  of  fragmenting  of  H+(H20)4 
clusters  which  are  initially  much  more  numerous  than  the  H+(H20)3  clusters. 
For  the  H+(H20)3  clusters,  fragmentation  is  not  a  50%  effect,  it  is  an 
effect  predicted  to  produce  an  increase  of  over  two  orders  of  magnitude. 
This  is  shown  quantitatively  in  Fig.  13,  which  the  integrated  effect  of 
both  fragmentation  and  agglomeration  on  H+(H20)3,  taking  into  account  the 
increase  brought  on  by  the  fragmenting  of  H*(H20)4  clusters.  (It  should  be 
noted  that  the  H+(H20)3  clusters  were  not  allowed  to  fragment  in  the 
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Figure  7.  The  Integrated  Effect  of  Agglomeration  on  Ion  Clusters  as  a 
Function  of  Axial  Position  for  the  30  Km  Test  Case. 


42 


INTEGRATED  LOSS  FROM  AGGLOMERATION 
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INTEGRATED  GRIN  AND  LOSS  FOR  H*(H20)3 


Figure  13.  The  Total  Calculated  Effect  of  Fragmentation  and  Agglomeration 
on  the  H+ (i^Olg  Number  Density  as  a  Function  of  Axial  Position  for  the  Two 
Test  Cases. 
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simulation,  since  no  H+(H20)2  clusters  were  carried  in  the  solution.  This 
process  would  not  have  a  substantial  effect  on  the  conclusion,  however, 
since  at  most  50*  or  so  of  the  H+(H20)3  clusters  would  fragment.  This 
would  still  leave  their  density  almost  two  orders  of  magnitude  above  what 
it  would  be  without  the  influence  of  fragmenting  H+(H20)4  clusters.) 

Figures  14  and  15  show  the  analogous  plots  for  H+(K20)4  and  H+(H20)5. 
In  the  former  case  the  net  effect  of  fragmentation  is  about  a  10-20* 
increase  in  number  density  and  in  the  latter  case  it  is  a  decrease 
consistent  with  Figs.  9  and  10.  This  is  because  the  density  of  H+(H20)6 
clusters  is  quite  low,  so  there  is  no  significant  source  term  for  H+(H20)5. 

6.4  Discussion 

Any  conclusion  from  these  calculations  is,  of  course,  provisional  in 
the  absence  of  experimental  verification.  On  the  basis  of  these 
calculations,  however,  the  two  major  conclusions  are  evident: 

1)  Agglomeration  of  additional  water  molecules  onto  ionic  clusters 
does  not  seem  to  be  a  significant  problem  in  mass  spectrometric 
sampling  of  stratospheric  ion  clusters.  This  is  primarily  due 
to  the  low  concentration  of  H2O. 

2)  Fragmentation  of  clusters  can  be  a  significant  problem.  It  is 

expected  that  a  significant  fraction  of  the  incoming  clusters 
will  undergo  fragmenting  for  the  experimental  conditions 
outlined  here.  For  clusters  such  as  where  there  are 

many  initial  clusters  with  more  water  molecules,  the  source  due 
to  the  fragmentation  of  larger  clusters  can  overwhelm  the 
original  population. 

Although  resources  did  not  allow  further  investigation,  it  is  quite 
possible  that  modifications  in  the  experimental  design  could  provide  less 
severe  fragmentation.  Clearly,  a  lowering  of  the  electric  field  strength 
would  produce  this  result,  but  it  would  be  at  the  cost  of  signal  strength. 
Another  solution  is  suggested  by  looking  at  Figs.  6  and  11.  It  can  be  seen 
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|  Figure  15.  The  Total  Calculated  Effect  of  Fragmentation  and  Agglomeration 
£  on  the  H+(H20>5  Number  Density  as  a  Function  of  Axial  Position  for  the  Two 
.*  Test  Cases . 


that  fragmentation  is  predicted  to  have  its  largest  effect  relatively  near 
the  orifice  where  the  collision  frequency  is  largest.  Figure  6  predicts, 
for  instance,  that  by  about  0.7  cm  from  the  orifice  the  fragmentation  rate 
is  an  order  of  magnitude  down  from  its  peak.  It  is  likely  that  potential 
grids  placed  beyond  this  point  would  produce  substantially  less 
fragmentation  while  maintaining  a  reasonable  signal  strength. 


52 


REFERENCES 


1.  Bird,  G.  A.,  Molecular  Gas  Dynamics,  Clarendon  Press,  Oxford,  1976. 

2.  Elgin,  J.  B. ,  "Monte  Carlo  Calculations  of  Mass  Spectrometer  Flow", 
Report  AFGL-TR-83-0057.  Air  Force  Geophysics  Laboratory,  February  1983. ADA128069 

3.  Shapiro,  A.H.,  The  Dynamics  and  Thermodynamics  of  Compressible  Fluid 
Flow,  the  Ronald  Press  Co.,  New  York,  1953,  pp.  73-105. 

4.  Sharafudinov,  R.G.  and  Skovorodko,  P.A.,  "Rotational  Level  Population 
Kinetics  in  Nitrogen  Free jets,"  Proceedings  of  the  12th  International 
Symposium  on  Rarefied  Gas  Dynamics.  AIAA  Press,  1980. 

5.  Bird,  G.A.,  "Breakdown  of  Continuum  Flow  in  Freejets  and  Rocket 

Plumes,"  Proceedings  of  the  12th  International  Symposium  on  Rarefied  Gas 
Dynamics,  AIAA  Press,  1980. 

6.  Liepmann,  H.W. ,  "Gaskinetics  and  Gasdynamics  of  Orifice  Flow,"  Journal 
of  Fluid  Mechanics,  Vol.  10,  No.  5,  1961,  pp.  65-79. 

7.  Smetana,  F.O.,  Sherrill,  W.A. ,  and  Schort,  D.R.,  "Measurements  of  the 

Discharge  Characteristics  of  Sharp-edged  and  Round-edged  Orifices  in  the 

Transition  Regime,"  Proceedings  of  the  5'th  Rarefied  Gas  Dynamics 
Symposium,  Vol.  2,  Academic  Press,  1967. 

8.  Labowski,  M.  Ryali,  S.,  and  Fenn,  J.B.,  "Flowfield  Calculations  in 

Nonequilibrium  Freejets  by  the  Method  of  Characteristics,"  Proceedings  of 
the  12th  International  Symposium  on  Rarefied  Gas  Dynamics,  AIAA  Press, 

1980. 

9.  Bird,  G.A.,  "Breakdown  of  Translational  and  Rotational  Equilibrium  in 
Gaseous  Expansions,"  AIAA  Journal ,  Vol.  8,  No.  11,  November  1970,  pp. 
1997-2003. 

10.  Bird,  G.  A.,  "Monte-Carlo  Simulation  in  an  Engineering  Context", 
Proceedings  of  the  12th  International  Symposium  on  Rarefied  Gas  Dynamics, 

Vol.  74,  Progress  in  Astronautics  and  Aeronautics,  AIAA,  New  York,  1981. 

11.  Vincenti,  Walter  G.,  and  Kruger,  Charles  H..  Jr.,  Introduction  to 
Physical  Gas  Dynamics,  John  Wiley  and  Sons,  1965,  pp.  348-356. 


12.  Borgnakke,  Claus,  and  Larsen,  Paul  S.,  "Statistical  Collision  Model 
for  Monte  Carlo  Simulation  of  Polyatomic  Gas  Mixture".  Journal  of 
Computational  Physics,  Vol .  18,  1975,  pp.  405-420. 

13.  Arnold,  F.,  Henschen,  G.,  and  Ferguson,  E.  E.,  "Mass  Spectrometric 
Measurements  of  Fractional  Ion  Abundances  in  the  Stratosphere-Positive 
Ions,"  Planet.  Space  Sci.  29(1981),  pp.  185-193. 

14.  Good,  A..  Durden,  D.  A.,  and  Kebarle,  P. ,  "Ion  Molecule  Reactions  in 
Pure  Nitrogen  and  Nitrogen  Containing  Traces  of  Water  at  Total  Pressures 
0.5-4  torr.  Kinetics  of  Clustering  Reactions  Forming  H+(H20)n,"  J.  Chem. 
Phys.  52(1970),  pp.  212-221. 

15.  Good,  A.,  Durden,  D.  A.,  and  Kebarle,  P.,  "Mechanism  and  Rate 
Constants  of  Ion-Molecule  Reactions  Leading  to  Formation  of  H+(H20)n  in 
Moist  Oxygen  and  Air,"  J.  Chem.  Phys.  52(1970),  pp.  222-229. 

16.  Lau,  Y.  K.,  Ikuta,  S.,  and  Kebarle,  P. ,  "Thermodynamics  and  Kinetics 
of  the  Gas-Phase  Reactions:  H-»0+{Ho0)n_..+H90s*Hq0+(H90)n,"  J.  Am.  Chem. 
SOC.  104(1982),  pp.  1463-1469.  3  2  "  1  2  3  2  " 

17.  See,  for  example,  Forst,  W.,  Theory  of  Unimolecular  Reactions 
(Academic  Press,  N.Y.,  1973). 

18.  Brady,  J.  W.  ,  Doll,  J.  D.  ,  and  Thompson,  D.  L.,  J.  Chem.  Phys. 
73(1980),  pp.  2767-2772. 

19.  Castleman,  A.  W.  ,  "Nucleation  and  Molecular  Clustering  About  Ions," 
Advances  in  Colloid  and  Interface  Science  H*(1979)«  PP-  73-128,  and 
references  therein. 

20.  Hsieh,  E.  T.  -Y.  and  Castleman,  A.  W. ,  "A  Reconsideration  of  the 
Theory  of  Capture  Cross  Sections  for  Ion/Molecule  Reactions  and  a  Total 
Energy  and  Angular  Momentum  Conserved  Average  Charge-Dipole  Interaction 
Theory,"  International  Journal  of  Mass  Spectrometry  and  Ion  Physics 
40(1981).  pp.  295-329. 

21.  Su,  T.  and  Bowers,  M.  T.,  J.  Chem.  Phys.  58(1973),  pp.  3027-3037.  See 
Ref.  20  for  criticism  and  corrections  concerning  this  approach. 

22.  Deguchi,  S.  and  Muhleman,  D.  0.,  "Mesospheric  Water  Vapor” .Journal  of 
Geophysical  Research,  87(C2),  Feb.  20,  1982,  pp.  1343-1346. 


23.  Kebarle,  P.,  Annual  Review  of  Physical  Chemistry,  28,  1977,  p.  445. 


Ballanthin,  J..  Private  Coaaunication. 

Handbook  of  Chemistry  and  Physics,  55th  Edition,  CRC  P’ess,  pp.  F58. 


